Models for diffusion and island growth in metal monolayers

We present a model which describes diffusion and island growth on FCC(001) metal substrates during molecular beam epitaxy. This model incorporates the essential characteristics of adatom hopping. To test it we performed Monte Carlo simulations for the growth of Cu on Cu(001) using semi-empirical potentials. We find that due to high edge mobility, islands are compact in the entire temperature range accessed by recent scanning tunneling microscopy experiments. The model predicts substantial dimer and trimer mobility and fast mobility of vacancies. Typeset using REVTEX 1 Recent experiments on thin film growth using molecular beam epitaxy have provided detailed information about growth kinetics and morphology. In particular the submonolayer regime has been studied extensively using both scanning tunneling microscopy (STM) and diffraction methods [1–8]. It was observed that for a variety of systems and a broad temperature range island nucleation is the dominant mechanism for crystal growth. A variety of island morphologies has been found, from fractal islands in Au on Ru(111) [2] to compact islands in Ni on Ni(001) [5]. Scaling properties of the island density and size as a function of deposition rate and coverage have also been studied [1–8]. These experiments stimulated much interest in theoretical studies of diffusion, island nucleation and morphology on surfaces using Monte Carlo (MC) simulations. In simulations of island growth during deposition atoms are deposited randomly at rate r and then hop according to a microscopic model [9–16]. Most studies use variations of the bond counting (BC) model in which the energy barrier for hopping into a vacant nearest neighbor site is given by EB = E0 + n · ∆E where E0 is the hopping barrier for an isolated atom, n is the number of nearest neighbors of the hopping atom in its initial site and ∆E is the bond breaking energy. The hopping rate is then determined by h = ν · exp(−EB/kBT ), where ν is the internal vibration frequency of the substrate, EB is the energy barrier, kB is the Boltzmann constant and T is the temperature. MC simulations show that this model exhibits both fractal islands, for large ∆E and more compact islands for small ∆E [12]. A related model was used to explore the island size distribution and its dependence on the critical island size [16]. The island morphology was also studied and it was shown that fractal islands are likely to appear for a broad temperature range on triangular substrates but not on square substrates, in agreement with STM results [15]. A different approach, introduced by Voter [17], was taken in Ref. [14] in which island growth of Cu on Cu(001) was simulated. In this study the hopping rates were determined by a complete set of energy barriers obtained using the atom embedding method of Finnis and Sinclair [18,19]. This approach can provide quantitative results but only limited understanding due to the large number of parameters. 2 In this letter we propose a simple model that has the simplicity of the bond counting model and at the same time it incorporates the essential physics of diffusion of adatoms on FCC(001) metal surfaces. In particular, the model shows high edge mobility which gives rise to compact island shapes, mobility of small islands such as dimers and trimers and high mobility of vacancies. The model is obtained from the analysis of microscopic energy barriers. We use a set of energy barriers for Cu on Cu(001) [20] obtained by the embedded atom method (EAM) [21]. This method uses semiempirical potentials and provides a good description of self diffusion on Cu(001) and similar surfaces [20]. Specifically, we use the EAM functions of Cu developed by Adams, Foiles, and Wolfer [22] which are fitted to a similar data base as the one employed by Foiles, Baskes, and Daw [23] except for the use of a more updated value of the single vacancy formation energy in the fitting [24]. The hopping energy barriers are calculated for all local environments in Fig. 1(a), where seven adjacent sites, i = 0, . . . , 6 are taken into account [25]. Each one of these sites can be either occupied (Si = 1) or vacant (Si = 0), giving rise to 2 7 = 128 barriers. A binary representation is used to assign indices to these barriers where E B, n = 0, . . . , 127 and n = ∑ 6 i=0 Si · 2 . The barrier height distribution is given in Fig. 2(a). We observe that this distribution features four peaks. When an atom on the surface hops into a vacant nearest neighbor site it has to cross the energy barrier between the initial and final sites. The top of this barrier is on the bridge site. Therefore, the hopping energy barrier is the difference between the energy on the bridge site and in the initial site. The occupancy of the 7 adjacent sites (Fig. 1(a)) affects both energies. We will now express the energy of the hopping atom in its initial site as: Ein = E 0 in −∆Ein · (S1 + S3 + S5)−∆E nnn in · (S0 + S2 + S4 + S6) (1) where Si = 1 if site i is occupied and 0 otherwise. The energy of an isolated atom is E 0 in while the reduction in its energy due the presence of an atom in a nearest (next nearest) neighbor site is given by ∆Ein (∆E nnn in ). Here we assume that the contributions of nearest neighbor and next nearest neighbor atoms to the energy are additive. This has been verified 3 by comparing the barriers obtained with the model and the EAM barriers, as described below. The energy of the hopping atom when it is on the bridge site is given by: Etop = E 0 top −∆Etop · (S1 + S2 + S5 + S6) (2) where E top is the energy of an isolated atom on top of a bridge site, while ∆Etop is the change in the energy due to the presence of an atom in one of the four sites adjacent to the bridge site. We do not include here n.n.n. type contributions since their effect is small. Therefore, the barrier EB = Etop − Ein for an atom to hop into an adjacent vacant site is given by: EB = E 0 B −∆Etop · (S2 + S6)+∆Ein · S3 − (∆Etop −∆Ein) · (S1 + S5) +∆E in · (S0 + S2 + S4 + S6) (3) where E B = E 0 top − E 0 in. To examine the formula we first found the parameters which best describe the 128 EAM barriers by minimizing the sum of squares ∑ 127 n=0[E n B(EAM) − E B(Eq. 3)] . The values obtained are E B = 0.494, ∆Ein = 0.265, ∆Etop = 0.268 and ∆E in = 0.024. We thus find that to within about 0.003 eV , ∆Ein = ∆Etop. Replacing both by ∆E = ∆Ein = ∆Etop we obtain a three parameter model (model I) in which the energy barrier is: EB = E 0 B +∆E · (S3 − S2 − S6) + ∆E nnn in · (S0 + S2 + S4 + S6). (4) The energy barrier distribution obtained from Eq. (4) is shown in Fig. 2(b). There is a good agreement in the location of the peaks with the EAM data, however the EAM peaks are significantly broader. One can further simplify the model by choosing ∆E in = 0. In the resulting model (model II), which has only two parameters (Fig. 1(b)), EB = E 0 B +∆E · (S3 − S2 − S6). (5) and the optimal values are found to be: E B = 0.526 and ∆E = 0.255. This model has the same complexity as the bond counting model (Fig. 1(c)), but it provides a far better description of diffusion on the Cu(001) substrate. 4 To examine our models we performed MC simulations of island growth for a range of deposition rates using both the EAM barriers and the two models [26]. We used a continuous time MC technique [14,27] in which moves are selected randomly from the list of all possible moves at the given time with the appropriate weights. The time is then advanced according to the inverse of the sum of all rates. The existence of four peaks in the spectrum of energy barriers indicates that there are four typical time scales of hopping. The two lowest peaks include very fast moves along island edges and motion of vacancies. The single atom move is in the third peak while moves in which atoms detach from islands are in the highest peak. From statistics collected during the simulations we find that most of the computer time is consumed by moves 6 (96) and 7 (112) which occur for atoms hopping along straight island edges. The reason is that for these moves the reverse move typically has the same low barrier. To make the simulation feasible we had to somewhat suppress these moves by artificially raising the barriers in the first and second peaks. This was done for all barriers lower than 0.4 eV according to: E B → E n B + α(0.4− E n B) where α = 0.7. Since the moves associated with these barriers are still orders of magnitude faster than for the two higher peaks we expect that this modification will have only small effect on the island morphology. We tried various values of α and found that up to α ≈ 0.8 this is indeed the case. In Fig. 3 we show the island morphologies obtained in MC simulations using both the EAM barriers and model I on a 200×200 lattice and compare them. The first comparison is at T = 270K, r = 0.0025 ML/s between the EAM barriers (Fig. 3(a)) and model I (Fig. 3(b)), and the second comparison at T = 230K, r = 0.01 ML/s between EAM barriers (Fig. 3(c)) and model I (Fig. 3(d)). In both cases there is very good agreement in the island morphology between the EAM barriers and model I. Model II also provides good qualitative agreement with EAM but island density is typically somewhat higher. We expect our model to apply also for other metal FCC(001) surfaces such as Ni(001) in which diffusion occurs through hopping rather than by the exchange mechanism. Preliminary EAM calculations for Ni(001) are consistent with our models. For surfaces in which the exchange mechanism is favorable, as there are indications for Al, Pt, Pd and Au modifications of the models are required. 5 However, we believe that the approach presented here should still provide a simple equation for the barriers as a function of the occupancy of the adjacent sites. Our models provide new and important insight about (a) the mobility of small islands (which has a strong effect on the island density and size distribution); (b) edge mobility (which determines the island shapes) and (c) vacancy mobility (which is the dominant diffusion mechanism in sputtering experiments). Small island mobility: in models I and II the diffusion coefficient for dimers (trimers) is approximately one half (one quarter) of the monomer diffusion coefficient with the same activation energy E B. This result is obtained by summing up all the possible moves of a dimer (trimer) and their rates. Therefore, at any temperature in which atoms are mobile, dimers and trimers are also mobile. This is radically different from the bond counting model where the barrier for dimer and trimer mobility is E B +∆E and thus typically turns on at higher temperatures than the single atom mobility. These predictions can be tested experimentally using field ion microscopy [28]. To examine the effect of small island mobility on the island density and size distribution we performed MC simulations of model II and the BC model. For both models we used the parameters E0 = 0.5eV , ∆E = 0.1eV and T = 250K and a range of deposition rates between 0.1 ML/s and 0.001 ML/s. We find that for this parameter range the island density for model II is about a half compared to the BC model and the peak in the island size distribution is shifted to the right [29]. Edge mobility: In our models the edge mobility is higher than the single atom mobility, since the barrier for hopping along a straight edge is E B−∆E compared to E 0 B for the single atom. This indicates, using an argument from Ref. [15], that for the entire temperature range in which islands appear they should be compact, unlike the case for the BC model [12]. This conclusion is confirmed by our MC simulations and is in agreement with recent STM results [5]. Note also, that in the BC model edge mobility is accompanied by detachment which has a similar barrier [12], while in our models there is no detachment in the relevant temperature range. Vacancy mobility: in models I and II vacancies are much more mobile than single atoms with 6 a barrier of the type E B − ∆E. Based on our MC simulations we predict that sputtering experiments in which single vacancies are created and then hop and nucleate should create much larger vacancy islands than the islands created in deposition experiments done under similar conditions. Previous simulations using a variation of the BC model lead to the opposite conclusion, namely that the vacancy islands are smaller and their density is higher compared to ordinary islands [11]. This is due to the fact that in the BC model vacancies have very low mobility and the energy barrier for their hopping is E B + 3∆E. In summary, we have proposed and studied two models that for the first time provide good quantitative description of diffusion and island nucleation on FCC(001) metal surfaces. Our models differ from the commonly used BC models on three aspects: dimer and trimer mobilities are comparable to the single atom mobility and the island density is thus reduced; the high edge mobility gives rise to compact island shapes in the entire temperature range of island growth; vacancy mobility is much higher than atom mobility. Since our models have only few parameters one needs a very small set of calculated energy barriers to determine them. This may open the way for a more effective use of first principle calculations of energy barriers on the surface as an input to the kinetic calculations. We would like to acknowledge support from the NSF under grants DMR-9217284 (O.B) and DMR-9119735 (G.V).